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ABSTRACT 

An analytic approximation to the mass function for gravitationally bound 
objects is presented. We base on the Zel'dovich approximation to extend 
the Press-Schechter formahsm to a nonspherical dynamical model. A simple 
extrapolation of that approximation suggests that the gravitational collapse 
along all three directions which eventually leads to the formation of real 
virialized objects - clumps occur in the regions where the lowest eigenvalue of 
the deformation tensor, A3, is positive. We derive the conditional probability 
of A3 > as a function of the linearly extrapolated density contrast 6 and 
the conditional probability distribution of S provided that A3 > 0. These two 
conditional probability distributions show that the most probable density of the 
bound regions (A3 > 0) is roughly 1.5 at the characteristic mass scale M*, and 
that the probability of A3 > is almost unity in the highly overdense regions 
{5 > 3a). Finally an analytic mass function of clumps is derived with a help of 
one simple ansatz which is employed to treat the multistream regime beyond 
the validity of the Zel'dovich approximation. The resulting mass function is 
renormalized by a factor of 12.5, which we justify with a sharp k-space filter by 
means of the modified Jedamzik analysis. Our mass function is shown to be 
different from the Press-Schechter one, having a lower peak and predicting more 
small-mass objects. 

Subject headings: cosmology: theory — large-scale structure of 

universe 
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1. INTRODUCTION 

Our present universe is observed to be quite clumpy with numerous galaxies, groups 
of galaxies, galaxy clusters, and etc. which span a large dynamic range in mass. The 
mass distribution function of these large scale structures is the crucial key to the nature of 
primordial density fluctuations from which the cosmic structures are believed to have arisen 
through gravitational growth, recollapse, and virialization ([Kolb fc Turner 1990|) . Since 
these gravitational processes are inherently nonlinear and sufficiently complicated, it is not 
an easy task to find the mass distribution function for bound objects analytically. Owing to 
its important role in cosmology, however, much effort has been made on determining even 
an approximate expression of the mass function (e.g., Peebles 1985| ; [Williams et al. 1991| ; 



Brainerd fc Villumsen 19"92| ; [Cavaliere fc Menci 1994| ; [Vergassola et al. 1994| ; Cavaliere, 



Menci, & Tozzi 1996). For recent review, see Monaco (1997). 

The pioneering attempt in this field has been ascribed to Press & Schechter (1974, 
hereafter PS) Q who proposed an analytic formalism for the mass function based on two 
simple assumptions: 1) the initial density field is Gaussian; 2) the gravitational collapse 
of mass elements is spherical and homogeneous. Along with these two assumptions, PS 
also postulated that the number densities of bound objects could be counted by filtering 
the initial linear density field. Although much criticism thereafter was poured upon the PS 
formalism about its unrealistic treatment of the collapse process and unclarified arguments 
including the notorious normalization factor of 2, the PS mass function has survived many 
numerical tests, showing good agreement with the results from N-body simulations (e.g., 
b]fstathiou et al. 1988| ; [Bond et al. 199T1 ; |Lacey fc Cole 1994|) . 



Motivated by somewhat unexpected success of the PS mass function, many authors 



^See also Doroshkevich (1967). 
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have tried to understand why it works so well in practice. Peacock & Heaven (1990) and 
Bond et al. (1991) have shown, by using the excursion set theory, that the fudge factor 
of 2 in the PS formalism which is directly related to the cloud-in- cloud problem can be 
justified with a sharp k-space filter. Jedamzik (1995) solved this cloud-in-cloud problem by 
means of the integral equation for the mass function. He insisted that the PS mass function 
should be altered even in the case of a sharp k-space filter. Yet Yano, Nagashima, & Gouda 
(1996) have argued that the sharp k-space filter recovers the PS mass function with the 
normalization factor of 2 even in the Jedamzik formalism if a mathematically consistent 
definition of isolated bound objects is used and the spatial correlations are neglected. 
They have also shown by introducing the two point correlation function into the Jedamzik 
formalism that the possible overlapping effect of density fluctuations which is responsible 
for the fragmentation and the coagulation of bound objects (see [Silk fc White 1978| ; [Lucchin 



19881 ; [Cavaliere fc Menci 1993|) can be neglected either on very small or on large mass scales. 



The PS approach has been also applied to nonspherical dynamical models. Monaco 
(1995) has suggested that the mass function should be treated as a Lagrangian quantity. 
Employing the Zel'dovich approximation as a proper Lagrangian dynamics, he computed 
the collapse epoch along the first principal axis, and showed that the shear shortens the 
collapse time and thus more high-mass structures are expected to form than the original PS 
mass function predicts. This effect of the shear explains dynamically the lowered density 
threshold {6c — 1.5) detected in several N-body experiments (e.g., pUfstathiou fc Rees 1988| ; 
Uarlberg fc Couchman 1989| ; |Klypin et al. 1995| ; |Bond fc Myers 1996|) . 



Shandarin fc Klypin (1984) have shown by N-body simulations that nonlinear clumps 
form from the Lagrangian regions where the smallest eigenvalue of the deformation tensor, 
A3, reaches a local maximum. Recently, Audit, Teyssier, fc Alimi (1997) have proposed some 
analytic prescriptions to compute the collapse time along the second and the third principal 
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axes, pointing out that Lagrangian dynamics is not valid after the first axis collapse but 
the formation of real virialized clumps must correspond to the third axis collapse. Their 
argument agrees with the N-body result obtained by Shandarin & Klypin (1984). In their 
analysis, the shear delays the third axis collapse rather than fastens it in contrast to its 
effect on the first axis collapse, which is in agreement with Peebles (1990). 

The normalization problem, however, has not been well addressed in these nonspherical 
approaches to the mass function. Monaco (1995) adopted the normahzation factor of 2 used 
in the PS formalism, while Audit et al. (1997) just assumed that the mass function could 
be normalized properly in any case. 

In this paper we study the eventual formation of clumps in a spatially fiat matter- 
dominated universe, with fragmentation and coagulation effects ignored. In § 2 we review 
the statistical treatment of the mass function, highlighting the PS formalism. In § 3 a 
nonspherical approach to the collapse condition based on the Zel'dovich approximation is 
described, and two useful conditional probability distributions relating the density field 
to the collapse condition are derived. In § 4 an ansatz is proposed to extend the validity 
of the Lagrangian dynamics to the third axis collapse. With a help of this ansatz an 
analytic approximation to the mass function for clumps is derived. In § 5 we justify the 
normalization factor 12.5 of the resulting mass function by using the Jedamzik integral 
equation. In § 6 the results are discussed and final conclusions are drawn. We relegate the 
detailed calculations and derivations to two Appendices. 



2. STATISTICAL DESCRIPTION OF MASS FUNCTIONS 

The mass function n{M) is defined such that n{M)dM is the comoving number 
density of gravitationally bound objects in the mass range (M, M + dM) . To compute this 
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statistics, it is assumed that the number densities of bound objects can be inferred from the 
hnearly extrapolated density contrast field, 6 = Sp/p (p: mean density). In other words, if 
a given region of the linear density field satisfies a specified criterion of collapse, then it is 
supposed to collapse and form a bound object. 

Let F{M) be the probability of finding a region satisfying a given collapse condition in 
the linear density field filtered at mass scale M, or equivalently the fraction of the volume 
occupied by the regions which will eventually collapse into bound objects with masses 
greater than M. Then we may write F{M) as follows: 

/oo 
piS) ■ CdS. (1) 
-oo 

Here p{6)dS is the probability that the smoothed density field at any given point will have 
a value in the range {6, 6 + d6), and C stands for the probability that the chosen point with 
density 6 will actually collapse. Once p{6) and C are determined, and then F(M) is found, 
the mass function n{M) can be easily obtained as 



dF 



dM ' 

where M/ p is nothing but the volume of a bound region with mass M. 



(2) 



The specific functional form of C is determined by the chosen dynamics to explain the 
collapse process, while p{6) depends on the property of the initial density field which is 
often assumed to be Gaussian in the standard cosmology (see [Bardeen et al. 1986[). The 



probability distribution of the Gaussian density field smoothed out by a window function 
W{R) of scale radius R is given by 



2a2(M) 

Here the mass variance cr^(M) is a function of scale mass M oc pR^ and estimated by 

d^k 



(3) 



^^M)= (4) 
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where Sk and Wk{R) are the Fourier components of the density S and the window function 
W{R) respectively. 

According to the top-hat spherical model adopted by PS, the bound objects form in the 
regions where the linearly extrapolated density contrast 6, growing with time, reaches its 
critical value 6c — 1.69 in a fiat universe (Peebles 1993|) . Therefore the regions with 6 > 6c 
when filtered at scale radius R correspond to the bound objects with masses greater than 
M{R) since it will have 6 = 6c when filtered at some larger scale. Thus in the PS formalism 
the collapse probability C in equation (1) is determined solely by the density field itself, 
and can be expressed by the following Heavyside step function: 

Cps = &i6-6c). (5) 

Using equations (1), (3), and (5), one obtains 

1 



F{M) 



1 



2na{M) J&c 





\ 1 


j exp 


[ 2a\M)\ 




\ 1 


/ exp 

he 


[ 2a^{M)\ 



-erfc 
2 



6c 



-V2a{M). 

where erfc(x) is the complementary error function. 



9(5 - 6c) d6 
d6 

(6) 



One obvious problem with the above analysis is that the integral of dF/dM over the 
whole range of mass does not give unity: 

/ ^dM = dF=-. (7) 
Jo dM Jo 2 ^ ' 

This normalization problem is originated from the fact that the PS formalism does 
not account for the underdense regions properly. Even for regions with 5 < (5c at a given 
filtering scale, there is still a nonzero probability that such regions will have 6 > 6c when 
filtered at some larger scale. But the PS formalism completely ignored those underdense 
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regions in estimating F{M), so half the mass initially present in the underdense regions was 
not taken care of. PS avoided this normalization problem simply by multiplying dF/dM by 
a factor of 2, and wrote the mass function in the form such that 



nps{M) 



'M 
ttM 



dF 


= 2^ 


da dF 


dM 


M 


dM da 



da 




dM 


a'^{M) 



exp 



2a2(M) 



As mentioned in § 1, the cooked up normalization factor of 2 in equation (8) has been 
shown to be correct in the case of a sharp k-space filter \Wk{R) = ©(vr/i? — k)], and various 
numerical tests have confirmed the PS mass function as a satisfactory approximation. 
Nevertheless it still leaves much to be desired: the physical meaning of the sharp k-space 
filter has yet to be understood; the gravitational collapse should be treated in more realistic 
models than the top-hat spherical one; the lowered density threshold {6c — 1.5) obtained in 
many numerical tests cannot be explained by this statistical argument, and so the PS mass 
function is degraded to a phenomenological device. 



3. NONSPHERICAL APPROACH TO COLLAPSE CONDITION 

Since PS derived their mass function on the basis of the top-hat spherical model in 
1974, the nonspherical nature of the gravitational collapse has been demonstrated by many 
authors (e.g., |Shandarin et al. 1995| ; Kuhlman, Melott, & Shandarin 1996). Especially 



the shear has been shown to play a very important role in the formation of the nonlinear 
structures (e.g., Peebles 1990| ; [Monaco 1995| ; [Audit et al. 1997| ). Therefore it is necessary to 



consider more realistic dynamical models to understand the collapse process and find the 
mass function. 

We choose the Zel'dovich approximation as a suitable Lagrangian dynamics to take 
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into account the nonspherical aspect of the gravitational collapse. However, instead of 
bringing the effect of the shear up to the surface, we try to retain the framework of the PS 
formalism, counting the number densities of bound objects from the filtered linear density 
field but with a different dynamical collapse probability C in equation (1). 



3.1. The Zel'dovich Approximation 

The Zel'dovich approximation (Zel'dovich 1970) asserts that the trajectory of a cosmic 
particle in the comoving coordinates can be expressed by the following simple formula: 

x = q-D+(t)Vvl>(q). (9) 

Here q and x are the Lagrangian (initial) and the Eulerian (final) coordinates of the particle 
respectively, \E'(q) is the perturbation potential which is a Gaussian random field, and 
D^(t) describes the growth of density fluctuations as a function of time. Throughout this 
paper, we focus on a spatially fiat matter-dominated universe with vanishing cosmological 
constant, in which case -D+(t) oc a{t) oc t"^^'^ [a{t): the cosmic expansion factor]. 

Applying a simple mass conservation relation pd^q = p(x)(i^x to the above formula (9) 
gives the following expression of the mass density: 

^^""^ ^ [1 - D4^)Ai(q)][l - D^[t)X,{q)][l - D+(t)A3(q)] ' ^^^^ 
where Ai, A2, A3 are the ordered eigenvalues (Ai > A2 > A3) of the deformation tensor, 

^ s^- 

Equation (10) shows that three random fields Ai(q), A2(q), A3(q) in the Lagrangian space 
are now the new dynamic quantities determining the collapse condition of given cosmic 
masses in the corresponding Eulerian space. Thus the mass function of bound objects can 
be built upon this Lagrangian dynamical theory (see also [Monaco 1995| ). 
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The actual dynamics for the formation of gravitationally bound objects is very complex. 
Even in the frame of the Zel'dovich approximation, the description of the gravitational 
collapse along all three directions is far from being simple and too cumbersome to use 
(Arnol'd, Shandarin, & Zel'dovich 1982). Here we employ rather a simplified dynamical 
model to approximate the collapse process and determine the collapse condition for the 
formation of clumps. 

Provided that at least one of the eigenvalues is positive at a given (Lagrangian) point, 
the denominator in equation (10) can become zero as -D+(t) increases with time, so the 
density p(x) will diverge, signaling collapse at the corresponding Eulerian point. If only 
the largest eigenvalue is positive (Ai > 0, A3 < A2 < 0) in a given region, then it collapses 
into a pancake. If two eigenvalues are positive (Ai > A2 > 0) while the third one is negative 
(A3 < 0), then a filament forms. The formation of a virialized bound object - a clump 
occurs only if all of three eigenvalues are positive, i.e. A3 > 0. So, in our dynamical model 
based on the Zel'dovich approximation, it is assumed that the lowest eigenvalue, A3, plays 
the most crucial role in determining the collapse condition for the formation of clumps. 
This assumption is in general agreement with Shandarin & Klypin (1984). 

The useful joint probability distribution of an ordered set (Ai,A2,A3) is derived by 
Doroshkevich (1970): 

P(Ai, A2, A3) = ( - ^ + ^) - - - ^12) 

where /i = Ai + A2 + A3, I2 = A1A2 + A2A3 + A3A1, and cr^ is the mass variance as defined in 
equation (4). Prom equation (12), one can see that the Zel'dovich approximation excludes 
both exactly spherical (Ai = A2 = A3) and exactly cylindrical (Ai = A2, A2 = A3, A3 = Ai) 
collapse. Both types of collapse have zero probability of occurring. (However, the points 
with Aj — Xj exist in generic fields on lines, that is on a set of measure zero in 3-dim., while 
the points with Ai = A2 = A3 do not exist at all.) 
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In order to obtain deeper qualitative understanding of the collapse in the Zel'dovich 
approximation, it may be also useful, in addition to this joint probability distribution (12), 
to have individual probability distribution of each eigenvalue (see Appendix A): 




The above individual probability distributions (13), (14), and (15) for the rescaled variable 
A/cr are plotted in Figure 1. Note that the distribution of A2(q) is Gaussian despite that 
A2(q) is not a Gaussian random field. 

According to equation (15), A3 > has a low probability of occurring, 0.08 (see 
Doroshkevich 1970| , or Appendix A). However, the small value of P(A3 > 0) = 0.08 does not 
indicate that only 8% of the whole regions will collapse into clumps. But rather it indicates 
that the probability of finding a bound region at filtering mass scale M is 0.08, provided 
that it is included in an isolated bound object with larger mass M' > M (see § 5). Here 
the isolated bound objects indicate the bound objects which have just collapsed at a given 
epoch. 

^Doroshkevich (1970) derived the probability distribution of Ai. But we found out a typo 
in his result. Except for the typo, equation (13) agrees with his result. 



In the following subsection, we derive the conditional probabilities of A3 > and 
5, reveal the correlated properties between them, and determine a nonspherical collapse 
probability C. 



3.2. Conditional Probabilities 

In the linear regime when D^{t) is still less than unity, equation (10) can be 
approximated by 

p~p[l + Z^+(Ai + A2 + A3)]. (16) 

Setting = 1 at the present epoch, the linearly extrapolated density contrast is now 
written as 

5 = ^ = Ai + A2 + A3. (17) 
P 

Let us choose (5, A2, A3) as a new set of variables. Then equation (12) can be reexpressed as 
a joint probability distribution of {6, A2, A3) such that 



+ ^(A2 + A3) (5 - A2 - A3) + ^AaAa 



x(5- 2X2 - A3) (A2 - A3) {5-X2- 2A3) . (18) 

Direct integration of the above joint distribution (18) over A2 gives the two point probability 
distribution of {S, A3): 

p{S,X3) = / ' p{5,X2,X3)dX2 

3\/5 A,, 2 nnx c .n.^2 o 2A f - 30A35 + 45Ai 

3^5 f 65^-30X38 + 45X1 



^i^'^'^K -)• ('^) 

where the upper limit and the lower limit of A2 are {5 — A3)/2 and A3 respectively due to 
the condition of Ai > A2 > A3. 
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With equation (19), we can investigate various correlated properties between 5-field 
and As-field. First of all, let us calculate the probability distribution of 6 confined in the 
regions with A3 > 0: 

p(5|A3>0) ^ PiSA3>0) _llp{6,Xs)dXs 



P(A3 > 0) P(A3 > 0) 

75^5 , / 
-0 exp 



25 / 52 



, — ■ exp , 

Here 6 stands for the Heavyside step function, and the condition Ai > A2 > A3 is used again 
to determine 6/3 for the upper limit of A3. Figure 2 compares the unconditional Gaussian 
distribution of the density field (3) with this conditional probability distribution (20) for 
the rescaled variable 6/a. It is shown that the maximum of p{S\X3 > 0) is reached when 
6 1.5a. That is, the linearly extrapolated density of the regions satisfying A3 > is most 
likely to be around 1.5cr. The average density contrast, < 6 >a3>0) can be also computed 
with equation (20): 



-H^)+erf(^)]}e(.). (20) 



< S >A3>o= / Sp{S\\s > 0)dS = -^-—{3^6 - 2) ^ 1.65a. (21) 

So in the regions with A3 > 0, the average density < 6 >a3>o is slightly higher than the 
most probable density, say S^^^q. We note that for cr = 1, S^^^q roughly coincides with 
the lowered density threshold 6c — 1.5 of the PS mass function, while < 6 >a3>o is close 
to the spherical threshold value 6c — 1.69. Setting a = 1 means filtering the density field 
on characteristic mass scale M^: [defined by (j(M*) = 1]. Thus the regions with A3 > for 
cr = 1 correspond to clumps with masses M > M*. In fact, as argued by Monaco (1995), 
it is unavoidable to limit our Lagrangian dynamical approach to the high-mass section 
(M > M^) since the Zel'dovich approximation is valid only in the single stream regions, 
while the multistream regions are rare for M > ( |Kofman et al. 1994| ). 
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Another conditional distribution worth deriving is P(A3 > 0\S), the probabihty that a 
given region with density 6 will have all positive eigenvalues: 

p{5, A3>0) 



P(A3 > 0\5) 



p{d) 



6 exp 



56^ 



The resulting conditional probability (22) for the rescaled variable 5/cr is plotted in Figure 
3. The probability of A3 > begins to exceed one-half when 6 ~ 1.5cr, and reaches unity 
when 6 ~ 3a. This implies that the collapse of highly overdense regions {6 ^ cr) will be 
always along all three directions (see also [Bernardeau 19941) . 

We take equation (22) as our nonspherical collapse probability C and proceed to derive 
the mass function of clumps analytically in the next section. 



4. AN ANALYTIC APPROXIMATION TO MASS FUNCTIONS 

As noted earlier, the Zel'dovich approximation as a first order Lagrangian theory works 
very well till the first orbit crossing (corresponding to the formation of pancakes) but breaks 
down afterwards in the multistream regime (phandarin fc Zel'dovich 1989|) . Therefore the 
rather restrictive collapse condition purely based on this Lagrangian formalism may not be 
fully satisfactory to describe the formation of clumps, especially low-mass objects. 

On the other hand, Shandarin & Klypin (1984) have shown by N-body simulations 
that the clumps form from the Lagrangian regions where the smallest eigenvalue A3 of the 
deformation tensor reaches a local maximum. Thus, one practical way to overcome the 
limited validity of the Zel'dovich approximation within the framework of our dynamical 
approach to mass functions is to parameterize the collapse condition by A3 > A3C, assuming 
that the critical value of A3C is a free parameter. Employing this simple ansatz to derive 
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n(M), we first calculate the following probability distribution with equations (3) and (19) 

s 

p{6,X3 > Age) _ Ixs,PiS,\3)d\3 



P(A3 > XsclS) 



p{6) 



p{6) 



{6 - SXsc) exp 



5(5-3A3c)' 



+ 2^ erf 



((5-3A 



3cy 



4a 



+ erf 



8a2 

(^-3A3,)v^ 
2a 



.e(<5-3A3c). (23) 



Comparison of equation (23) with equation (22) reveals that P(A3 > XsdS) is just 
horizontally shifted along (5-axis by 3A3c from P(A3 > 0|5) with its shape unchanged. 

Consequently, this ansatz is mathematically equivalent to parallel transformation of 
the density field itself by — 3A3c. Thus equation (1) is now expressed as follows: 



/oo 
p(5 + 3A3j-P(A3>0|5)d5 
-oo 



-oo 

1 



27ra Jo 



exp 



(5 + 3A 



3c 



+ - 



25 



2(t2 



75VT0, / 55\ 



(24) 



Here the volume fraction F{M) is normalized by a factor of 1/0.08 = 12.5 which we justify 
with a sharp k-space filter in § 5. This normalization factor is much larger than the factor 
of 2 in the PS formalism. However, this larger normalization factor can be explained by 
the larger amount of cloud-in- cloud occurrences in our dynamical formalism than in the PS 
formalism, as shown in § 5 where such amount is computed. In an ideal hierarchical model, 
all the masses are included in clumps. According to our dynamical model, only about 8% 
of all the masses are included in the clumps with the "largest" mass (the "largest mass" of 
bound objects in the universe is, in a practical sense, M ~ M^,). This is in rough agreement 
with the fraction of the galaxies in the Abell clusters (e.g., [Padmanabhan 1993| ). All the 
remaining masses are included in the clumps at smaller filtering mass scales. 
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Differentiating equation (24) witli respect to cr, we liave 

(5 + 3A3c)^r 75v^0 



dF _ d j 1 p / 



2a2 ^ I 80F(T 
'(5A/iO 



(5 exp 



55^ 
8(t2 

(5a/T0 



25V10A3C/5A 



20Fa2 



^3c 



H erf + erf 

4 [ V 4a ^ ^ 2a 

±)exp(-§)erfc(^) 



ci5 



, 25V15A3C 
H — ^ ^ , exp 



15A2 X ^ /v^A 
erfc 



3c 



125V5A^ 



3c 



exp 



9A 



3c 



(25) 



8v^a2 "^"''V 4^2 y--v 2a ^ 67ra3 "^""^ V 2a2 

Figure 4 sliows tlie generic beliavior of this differential volume fraction (25) as A3C changes. 
Since OF /da is directly proportional to n(M), one can conclude from Figure 4 that as A3C 
increases, the number densities of small-mass clumps (large a) increase while the large 
masses (small a) are reduced and the peak is lowered. 



For simple power law spectra oc fc", the mass variance becomes 



-(n+3)/3 



(26) 



So in this case, the mass function can be expressed explicitly in terms of M: 
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dF 
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da dF 
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dM 
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25v^A3c/n + 3\ p /M\("+3)/6 



2A 



5A2^ / M \ 



6 JM^\MJ 
1 
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X 



exp 



5A2^ 1^ M \("+3)/3 



erfc 



, \/6 

H — — exp 



15A2^ / M \ ("+3)/3- 



4 VM, 
5A3e /M\('^+3)/6 



erfc 



3V2^VM, 



exp 



9AL 



2 \MJ 
VSXsc ( M\ ("+3)/6 
2 \mI 

M \ ("+3)/3' 



V2X 



3c 



M \ ("+3)/6' 



(27) 



We display the resulting mass function for A3C — 0.37 in Figure 5. The value of 0.37 
for A3C is chosen to make our results for the high-mass section fit well with the 6c — 1-5 PS 



mass function which has been tested to be a good approximation (see [Monaco 1995|) . The 
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original PS mass function with 5^ — 1-69 is also shown for comparison. For every power 
index n from —2 to 1, the mass function (27) is characterized by the following properties: 

(1) In the high-mass section (M/M* > 1), it fits quite well with the Sc — 1.5 PS mass 
function. 

(2) Its peak is lower than that of the PS one, which agrees with N-body results (e.g., 
see [Efstathiou et al. 1988|) . 

(3) It has approximately the same slope as the PS mass function but predicts more 
structures in the low- mass section (M/M^, < 1). 

5. NORMALIZATION 

Up to now, following the PS-like approach, we assumed that equation (2) is correct. In 
other words, the probability of finding a region with A3 > X^c at filtering mass scale M is 
assumed to be proportional to the fraction of the volume occupied by the regions which will 
eventually collapse into bound objects with masses > M. 

Strictly speaking, however, equation (2) is not quite correct since the resulting mass 
function has to be always renormalized. Even in the regions with A3 < A3C at the filtering 
mass scale M, there is still a nonzero probability of A3 = A3C when the density field is filtered 
at some larger scale M'{> M). But this marginal probability is ignored by the PS-like 
approach, which has resulted in a large normalization factor 12.5 of our mass function. 
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Jedamzik (1995) suggested a generalization of equation (2) Q 



dF 



dM 



d 



dM 



oo 

dM'n(M')—P(M, M'] 
p 



(28) 



Here P{M, M') is the conditional probability of finding a bound region (A3 > Asc) at 
filtering mass scale M, provided that it is included in an isolated bound object (A3 = Asc) 
with mass M'(> M). The isolated bound objects at a given epoch are those which have 
just collapsed. Thus, in our formalism, the isolated bound objects correspond to the regions 
with A3 = A3C at a given filtering mass scale. 

We find that the conditional probability P{M, M') for the case of a sharp k-space filter 
is given by (see Appendix B): 



P(M, M') = 0.08e(M' - M). 



(29) 



Equation (29) reveals that P(A3 > 0) = 0.08 results in P(M, M') = 0.08 (M < M'). So, in 
our formalism the probability of finding a bound region (A3 > A3c) of mass scale M included 
in an isolated bound region (A3 = A3c) with mass greater M is only 0.08. And this is directly 
related to our normalization factor of 1/0.08 = 12.5. Whereas in the PS formalism the 
probability P{M, M') is 0.5, which is again directly related to the PS normalization factor 
of 1/0.5 = 2. 



Now, with equation (28) and (29), we have 



dF 



dM 



/■oo Ml M 

0.08 / dM'n(M')—6D(M' - M) = 0.08— n(M), 
Jo p p 



(30) 



where 5d stands for the Dirac delta function. More explicitly. 



n(M) = 12.5 



M 



dF 



dM 



(31) 



■^In the original analysis based on the top-hat spherical model, Jedamzik (1995) did not 
use a mathematically correct definition of isolated bound objects. 
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which is exactly the same formula as equation (2) with the normalization factor of 12.5 
included explicitly. 

Thus equation (31) justifies the normalization factor 12.5 of our mass function in the 
case of a sharp k-space filter. 

6. DISCUSSION AND CONCLUSIONS 

We have derived an analytic approximation of the mass distribution function for 
clumps. The underlying dynamics has been described by the Zel'dovich approximation 
which treats the nonspherical gravitational collapse. Similar to Shandarin & Klypin (1984) 
and Audit et al. (1997), we have assumed that the clumps would be formed by the mass 
elements which have experienced gravitational collapse along all three directions. 

We have given a somewhat different interpretation to the PS analysis by rccxpressing 
the fraction of the volume occupied by the bound regions in terms of two probabilities: 
the probability of the Gaussian density distribution and the collapse probability of the 
given dense regions. The Zel'dovich approximation has led us to determine a nonspherical 
collapse probability which is different from the PS one, relating the density field to the 
positive lowest eigenvalue of the deformation tensor. 

We have shown that the collapse probability reaches almost unity when S > Sa, which 
indicates that the highly over dense regions will always collapse along all three directions. 
In addition, we have found the density distribution of the regions which meet the collapse 
condition based on the Zel'dovich approximation. This distribution has shown that the 
most probable density contrast of such regions is around 1.5 at the characteristic mass scale 

We have proposed a simple ansatz in order to treat the multistream regions where 
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the Lagrangian dynamics is not applicable. This ansatz has enabled us to derive an 
analytic mass function characterized by one free parameter, X^^- The best approximate 
value of this parameter has been chosen to be Asc — 0.37. We admit that there is no 
background dynamical theory for determining directly the value of this free parameter 
which thus has to be found phenomenologically. However, the following arguments may 
give a hint to understand why this parameter has this value. A simple extrapolation of 
equation (10) into the multistream regions suggests that only the mass elements with 
A3 > Asc = 1 collapse along all three directions by the present epoch of = 1. However, 
the collapse along the first two directions increases the density, which therefore speeds 
up the collapse along the third directions. This roughly agrees with the conclusion of 
Audit et al. (1997). Using equation (24) from Audit et al. (1997) with their choice of the 
parameters {e = l,a = 0.8,5c = 1.69, ctc = 0.74), one can easily obtain for the collapse 
epoch Qc = 1/(0. 8A3 + 0.325) that is always earlier than the prediction of the Zel'dovich 
approximation (ac = I/A3) provided that A3 > 0. 

For power law spectra, it has been shown that our resulting mass function with 
Asc — 0.37 is in good agreement with the 6c — 1.5 PS mass function in the high-mass section, 
but has lower peak and predicts more small-mass structures, which are in agreement with 
what has been detected in N-body simulations (e.g., [Efstathiou et al. 1988| ; [Feacock~S 



Heavens 1990| ). However, it should be noted that the prediction concerning the small-mass 



structures is least reliable not only in any PS-like approach, but also in our dynamical 
approach to the mass function since the validity domain of the Zel'dovich approximation is 
limited to the high-mass section as outlined in § 3.2. 

Like the other PS-like formalisms, a normalization factor for the mass function has 
been introduced, which in our case is 12.5. We have justified the normalization factor with 
a sharp k-space filter by using the Jedamzik integral equation, showing that this rather 



- 21 - 



large normalization factor is due to the low probability of finding a bound region (A3 > X^c) 
at filtering mass scale M included in an isolated bound region (A3 = Asc) with larger mass 
M'. But the physical meaning of the sharp k-space filter has yet to be fully understood. 

We postpone the numerical testing of our mass function to the following paper. 

We are grateful to Lev Kofman, Paolo Catelan, and the referee for useful discussions 
and helpful comments. This work has been done under the support of NASA grant NAG 
5-4039 and EPSCoR 1996 grant. 
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A. APPENDIX 

In 1970 Doroshkevich found the joint probability distribution, p(Ai, A2, A3) of an 
ordered set of eigenvalues [equation (12) in § 3.1], corresponding to a Gaussian potential. 
In this appendix, we sketch the derivation of p(Ai), p{X2), and ^(As) [equation (13), (14) 
and (15) in § 3.1], and investigate their statistical properties. 



The two point probability distributions, p(Ai, A2),p(A2, A3) and p(Ai, A3) can be easily 
obtained from the direct integration of p{Xi, X2, A3) such that 



p(Ai, A2 



/A2 
p(Ai, A2, A3)(iA3 
-00 

(Ai - A2)(3Ai - A2) exp 



-00 
1125 

64"\/57r(7' 



3A^ ^ 3A1A2 



9Ai 



37r(T 



12 
X exp 



(Ai - A2) 
2 



8 + — (3Ai- A2)(3A2- Ai; 



45Af ^ I5A1A2 



16(t2 



8a2 



45A| 
16(t2 



erfc 



V3 

4(7 



(Ai - 3A2 



2(72 



.e(Ai-A2), (AI) 



P(A2, A3) 



A2 



p(Ai, A2, A3)(ZAi 



1125 



64V57r(j4 



(A2 - A3) (A2 - 3A3) exp 



3A3 3A2A3 



9A2 



X exp 



3 



8+^(A2-3A3)(A3-3A2) 



45A? I5A2A3 45A2 



16(72 



+ 



8(72 



16(72 



erfc 



^(3A2-A3) 



2(72 



.e(A2-A3), (A2) 



p(Ai, A3) 



Ai 



p{\l, A2, X3)dX2 



1125 



64\/57r(7^ 



(Ai - A3)(3Ai - A3) exp 



3A2 3A1A3 



9Ai 



a" 



2(72 



+(Ai - A3)(Ai - 3A3) exp 
3 



3A3 3A1A3 



9A| 
2(t2 



/37r(7,, , , 
+^^(Ai-A3) 



(3Ai-A3)(Ai-3A3)-8 
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X exp 



X < erfc 



16(72 
4(7 



8(72 



A3) 



45A| 
16(t2 



erfc 



^(3A3-A0 



}|e(Ai-A3) 



(A3) 



In order to derive p{\i), p{^2), and p{\3), we have to integrate the above two point 
distributions, which involve complex error function terms as one can see. We find the 
following recursion formula which is useful in integrating such complex terms. 



n 



2o? 
1 

h 



+ 



a' 



- yr-'exp(-aV)erf(bt)di 
exp(-a2t2)erf(bt) 
= /r-^exp[-(a' + &')i2]di. 

TT J 



(A4) 



With the above recursion formula, it is straightforward to derive the individual distributions. 
The results are shown in § 3.1 [equation (13), (14), and (15)]. 

Now the probability that each eigenvalue is positive as well as the mean and the 
variance of each eigenvalue can be computed with the above results: 



P(Ai > 0) 
A"i = 



23 
25' 



P{^2 > 0) = ^, P(A3 > 0) = ^, 



= (T, 



137r - 27 



307r 



A2 = 0, A3 ^ 

2 

(T\ = (T 



x/IOtt 



ISn - 27 



307r 



-cr 



(A5) 
(A6) 
(A7) 



which are all in agreement with Doroshkevich (1970). 



B. APPENDIX 



In the framework of our formalism, the conditional probability P{M, M') is written as 

P(M, M') = P(A3 > A3e| a;, = A3e) = ^(^3 > A3. A^ ^ A3c) ^ ^^^^ 
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where A3 and A3 are the lowest eigenvalue of the deformation tensor at the same point but 
at two different filtering mass scale M and M' respectively. 

In order to derive the probability P(A3 > Aac, A3 = Asc), we start with the multivariate 
Gaussian joint probability distribution (see Poroshkevich 1970| ; Pardeen et al. 1986|) for the 



six independent elements of the deformation tensor at two different mass scales: 

Pj{yi,---,y6)dyi---dy'(i = —^^^^^j===dyi---dy'Q, (B2) 

(27r)6y^det(V) 

Q = ^y'-V-y (B3) 

Here V is the covariance matrix, while {yi}^^i and {y'i}^^i are the six independent elements 
of the deformation tensor (defined by equation (11) in § 3.1; yi = dn, y2 = ^22, I/3 = c^ss, 
2/4 = di2, y5 = d23, ye = rfsi) at mass scale M and M' respectively. 

In the case of a sharp k-space filter, the mutual correlations between {yi}^^i and 
{y[}Li are 

<y^>=^, <y.y,>=^, <y^y>="\. m 
<v'^>=% <v\v]>=% <v^]>=% (B5) 

for i, j(7^ I) = 1, 2, 3, and 

2 2 
<y->=^, <y^yj>=0, <y,y[>=^, (B6) 

<yf>=^, <i/:i/;>=o, <i/.y;>=o, (b?) 

for i) = 4, 5, 6. Here afj and crfj/ are the mass variance of the density field filtered at 

the mass scale M and M' respectively. 

Through equation (B2) to (B7), along with the similarity transformation of the 
deformation tensor into its principal axes, we find the following joint probability distribution 
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of the three eigenvalues {Xi}f=n {Xi}i=i of the deformation tensor at filtering mass scale M, 
M' (M < M') respectively. 



Pj(Xi, ■ ■ ■ , X'^)dXi ■ ■ ■ dX's 



Pji(Ai, As, A3)dAidA2dA3 
xpj2{X[,X',,X',)dX[dX',dX', 



5^ -S^ 



Pji 

PJ2 



53-33 
2%V|^V5 



exp 



3/i2 15/2 



4 



2ai 



(Ai-A2)(A2-A3)(A3-Ai), 



exp 



3/f 15/^ 



+ 



M' 



(a;-a'2)(a'2-a'3)(a^-a;), 



(B8) 

(B9) 
(BIO) 



where 



a, = a,-a:, 

7i = Ai + A2 + A3, 72 = A1A2 + A2A3 + A3A1, 

I[ = X[ + A2 + A3, I2 = X[X2 + A2A3 + AgA^. 



(Bll) 
(B12) 
(B13) 



Note the similarity between pji and pj2. In fact the above equations hold good only in the 
case of a sharp k-space filter. 

The integration of pj over Ai, A2, A^, and A2 gives us the joint probability density 
distribution, p(A3, A3): 

p(A3, A'3)dA3dA'3 = p(A3)dA3p(A^)dA'3. (B14) 

Here the probability density distributions of p( A3) and ^(Ag) have the same form as p(A3) 
[equation (15)] except for the value of the variance. 

Finally we derive the conditional probability P{M, M'): 

P{Xs > A3e, A'3 = A3c) P(A'3 = A3e)c^A'3 J^^ dAscPiA^,) 



P{M, M') 



PiX's = A3c) 

dAs,p{As,) = 



p{X's = Xs,)dX', 
c^A3p(A3) = 0.08e(M' - M), 



(B15) 



where A3C is A3 — A 
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Fig. 1. — The individual probability distributions of three eigenvalues of the deformation 
tensor : The sohd lines shows the analytic results obtained in this paper for the rescaled 
variable A/u. The numerical results from the Monte Carlo simulation are also plotted as the 
dotted lines. 
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Fig. 2. — The probability distribution of the rescaled density field, {5 /a). The solid line 
represents the rescaled density distribution under the condition of A3 > 0, while the dashed 
line shows the unconditional Gaussian distribution. The vertical dotted line indicates the 
position of 5/(7 = 1.5 
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Fig. 3. — The conditional probability of A3 > as a function of the rescaled density S/a. 
The vertical dotted hne indicates the position S/a — 1.5. 
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Fig. 4. — The differential volume fraction for Asc — 0.1,0.4,0.7 and 1.0 (dot-dashed, soHd, 
dashed and long-dashed lines). The dotted line is the standard {6c — 1.69) PS differential 
volume fraction. 
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Fig. 5. — The mass function for the power index n = +1, 0, — 1 and —2. The sohd hne shows 
our analytic results with X^c — 0.37, while the dashed and dotted lines represent the PS 
mass function with Sc — 1.5 and 1.69 respectively. 



